# include <stdio.h>
# include <stdlib.h>
# include <unistd.h>
# include <time.h>
# include <math.h>
# include <gmp.h>

int main(void) {
    int i, k;
    mpf_t pi, a0, a1, y0, y1;
    mpf_t first_term, second_term, third_term, fourth_term;
    mpf_t aux;
        FILE *fp;
        fp = fopen("boweinG.txt","w+");

    mpf_set_default_prec (1000000);

    mpf_init(pi);
    mpf_init(a0);
    mpf_init(a1);
    mpf_init(y0);
    mpf_init(y1);
    mpf_init(aux);

    mpf_init(first_term);
    mpf_init(second_term);
    mpf_init(third_term);
    mpf_init(fourth_term);

    //a[0] = 6 - 4*sqrt(2);
    mpf_sqrt_ui(first_term, 2);
    mpf_mul_ui(first_term, first_term, 4);
    mpf_ui_sub (first_term, 6, first_term);
    mpf_set(a0, first_term);

    //y[0] = sqrt(2) - 1;
    mpf_sqrt_ui(first_term, 2);
    mpf_sub_ui (first_term, first_term, 1);
    mpf_set(y0, first_term);

    for(k=0;k<11;k++) {
        //y[k+1] = (1 - pow((1 - pow(y[k],4)),0.25))/(1 + pow((1 - pow(y[k],4)),0.25));
        mpf_pow_ui(first_term, y0, 4);
        mpf_ui_sub(first_term, 1, first_term);
        mpf_sqrt(first_term, first_term);
        mpf_sqrt(first_term, first_term);
        mpf_ui_sub(first_term, 1, first_term);
        mpf_pow_ui(second_term, y0, 4);
        mpf_ui_sub(second_term, 1, second_term);
        mpf_sqrt(second_term, second_term);
        mpf_sqrt(second_term, second_term);
        mpf_add_ui(second_term, second_term, 1);
        mpf_div(first_term, first_term, second_term);
        mpf_set(y1, first_term);

                   //first_term              //second_term    //third-term
        //a[k+1] = a[k]*pow((1+y[k+1]),4) - (pow(2,(2*k+3)) * y[k+1]*(1+y[k+1]+pow(y[k+1],2)));
        mpf_add_ui(first_term, y1, 1);
        mpf_pow_ui(first_term, first_term, 4);
        mpf_mul(first_term, a0, first_term);

        mpf_mul_ui(second_term, y1, (long) pow(2, (2*k)+3));

        mpf_pow_ui(aux, y1, 2);
        mpf_add(aux, y1, aux);
        mpf_add_ui(aux, aux, 1);
        mpf_mul(second_term, second_term, aux);
        mpf_sub(first_term, first_term, second_term);
        mpf_set(a1, first_term);

        //pi = 1/a[k+1];
        mpf_ui_div(first_term, 1, a1);
        mpf_set(pi, first_term);

        mpf_set(y0, y1);
        mpf_set(a0, a1);

        gmp_printf("%.50Ff\n", pi);
    }
mpf_out_str(fp, 10,0, pi);

    return 0;
}
